simulation_brooklyn<-function(model,n,smooth=5,bwsel="CCT"){

## fitted model from the Brooklyn Tech data
if (model==1) {
  c.p = c(0.56, 0.18, -0.11, 0.19)
  c.m = c(0.70, 0.43, 0.63, 0.38)
  u.rnd.mnus  = function(n) { rnorm(n)*0.2}
  u.rnd.plus = function(n) { rnorm(n)*0.2 }
} else
if (model==2) {
  c.p = c(0.61, -0.02, 0.06, 0.17)
  c.m = c(0.61, -0.02, 0.06, 0.17)
  u.rnd.mnus  = function(n) { rnorm(n)*0.2 }
  u.rnd.plus = function(n) { rnorm(n)*0.2 }
} else
if (model==3) {
  c.p = c(0.61, -0.02, 0.06, 0.17)
  c.m = c(0.61, -0.02, 0.06, 0.17)
  u.rnd.mnus  = function(n) { rnorm(n)*0.2 }
  u.rnd.plus = function(n) { rnorm(n)*0.05 }
}

  X = 2*rbeta(n,2,2)-1 
  u.l  = u.rnd.mnus(n) 
  u.r  = u.rnd.plus(n) 
  y.l =c.m[1] + c.m[2]*X + c.m[3]*X^2 + c.m[4]*X^3   + u.l
  y.r =  c.p[1] + c.p[2]*X + c.p[3]*X^2 + c.p[4]*X^3   +  u.r  
  y = c(y.l[X<0],y.r[X>0])
  x = c(X[X<0],X[X>0])

# MRD test
mrdbw.out = rdbwselect(y,x, bwselect = bwsel)
bw_mrd=mrdbw.out$bws[1,1]*n^(1/5-1/smooth)
if (is.finite(bw_mrd)==1){
results<-MRD_sharp(y,x,bw_mrd,Xeval=0)
MRDbw=bw_mrd
MRD=results$estimate
MRDstat=results$tstat
# DRD test
ygrid=quantile(y,seq(0.001,0.999,length.out=999))
results<-DRDtest(y,x,ygrid,h=bw_mrd,bwsel=bwsel,undersmooth=smooth)
DRDbw=results$bw
DRDstat=max(abs(results$ksstat))
DRDbw2=results$bw_iter
DRDstat2=max(abs(results$ksstat_iter))
# confidence region
set1<-ygrid[results$ksstat_iter<1.2239]
set2<-ygrid[results$ksstat_iter< -1.2239]
if (model==1){
  rightsettrue=NA
  leftsettrue=ygrid
} else
if (model==2) {
  rightsettrue=NA
  leftsettrue=NA
} else
if (model==3) {
  rightsettrue=ygrid[ygrid< 0.61]
  leftsettrue=ygrid[ygrid> 0.61]
}
} 
else {
MRDbw=NA
MRD=NA
MRDstat=NA
DRDbw=NA
DRDstat=NA
DRDbw2=NA
DRDstat2=NA
}
list(MRDbw=MRDbw,MRD=MRD,MRDstat=MRDstat,DRDstat=DRDstat,DRDstat2=DRDstat2,DRDbw=DRDbw,DRDbw2=DRDbw2,set1=set1, set2=set2,leftsettrue=leftsettrue,rightsettrue=rightsettrue)
}

